In [1]:
%matplotlib inline
from functools import cache
import numpy as np
import matplotlib.pyplot as plt
import itertools
import math
from decimal import Decimal, getcontext
In [2]:
def make_pqzo(p1, p2, p3):
    @cache
    def p(i, k):
        if k < 0:
            raise ValueError("p called with k < 0")
        if k == 0:
            return [p1, p2, p3][i - 1]
        return p(i, k - 1) * (p(i, k - 1) + 2 * p(i % 3 + 1, k - 1))
    
    @cache
    def q_odd(i, k):
        return q(i, k - 1) * (q(i, k - 1) + 2 * q(i % 3 + 1, k - 1))
    
    @cache
    def z(i, k):
        return p(i, k) * (p(i, k) + 2 * p((i+1) % 3 + 1, k))
    
    @cache
    def q(i, k):
        if k < 0:
            raise ValueError("q called with k < 0")
        if k == 0:
            return p(i, 0) * (p(i, 0) + 2 * p((i + 1) % 3 + 1, 0))
        return (
            q_odd(i, k) * (z(i, k) + z(i % 3 + 1, k))
            + q_odd(i % 3 + 1, k) * z(i, k)
        )
    
    return (p, q, z, q_odd)
In [3]:
def plot_pq(p, q, z, o, k, **kwargs):
    plot_pq_range(p, q, z, o, 0, k, **kwargs)

def plot_pq_range(p, q, z, o, k_low, k_high, simplex=True, simplex2=True, basic=True, log_scale=False):
    if simplex:
        fig, ax = plt.subplots()
        X = np.array([1, -1, 0])
        Y = np.array([0,  0, 1])
        for fun, name in [[p, 'p'], [q, 'q'], [o, 'o']]:
            x_axis = []
            y_axis = []
            for k in range(k_low if fun != o else max(k_low, 1), k_high):
                point = np.array([fun(1, k), fun(2, k), fun(3, k)])
                x_axis.append(point @ X)
                y_axis.append(point @ Y)
            ax.plot(x_axis, y_axis, '->', label=name)
            #ax.plot(x_axis, y_axis, '-4', label=name)
        ax.set_xlim(-1.1, 1.1)
        ax.set_ylim(-.1, 1.1)
        ax.legend()
        fig.show(warn=False)
    
    if simplex2:
        fig, ax = plt.subplots()
        X = np.array([1, -1, 0])
        Y = np.array([0,  0, 1])
        
        # Plot p
        x_axis = []
        y_axis = []
        for k in range(k_low, k_high):
            point = np.array([p(1, k), p(2, k), p(3, k)])
            x_axis.append(point @ X)
            y_axis.append(point @ Y)
        ax.plot(x_axis, y_axis, '->', label="p")
        
        # Plot q
        x_axis = []
        y_axis = []
        for k in range(k_low, k_high):
            point = np.array([q(1, k), q(2, k), q(3, k)])
            x_axis.append(point @ X)
            y_axis.append(point @ Y)
            if k + 1 < k_high:
                p2 = np.array([o(1, k+1), o(2, k+1), o(3, k+1)])
                x_axis.append(p2 @ X)
                y_axis.append(p2 @ Y)
        ax.plot(x_axis, y_axis, '->', label="q")
        
        ax.set_xlim(-1.1, 1.1)
        ax.set_ylim(-.1, 1.1)
        ax.legend()
        fig.show(warn=False)
    
    if basic:
        fig, ax = plt.subplots()
        if log_scale:
            ax.set_yscale("log", base=2)
        colors = []
        for high, low in [(1, 0), (.6, .2)]:
            colors.extend([(high, low, low), (low, high, low), (low, low, high)])
        color_index = 0
        for i in [1, 2, 3]:
            data = []
            x_vals = []
            for k_val in range(k_low, k_high):
                x_vals.append(k_val)
                data.append(p(i, k_val))
            ax.plot(x_vals, data, label=f"$p_{i}$", color=colors[color_index])
            color_index += 1
        for i in [1, 2, 3]:
            x_axis = []
            y_axis = []
            for k in range(k_low, k_high):
                x_axis.append(k)
                x_axis.append(k + .5)
                y_axis.append(q(i, k))
                y_axis.append(o(i, k + 1))
            ax.plot(x_axis, y_axis, label=f"$q_{i}$", color=colors[color_index])
            color_index += 1

        ax.legend()
        # ax.set_title
        # fig.savefig
        fig.show(warn=False)
In [4]:
getcontext().prec = 2000
for i in range(1, 32, 2):
    a = Decimal(i) / Decimal(100)
    plot_pq(*make_pqzo(a, 1 - 2*a, a), 100)

<ipython-input-3-3d26527cd5ca>:56: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots()
In [5]:
getcontext().prec = 2000
epsilon = Decimal(i) / Decimal(2**10)
c = epsilon ** 2
a = epsilon * 2
b = 1 - a - c
plot_pq(*make_pqzo(a, b, c), 30) # note p1 high in range 10–30
In [6]:
getcontext().prec = 10000
epsilon = Decimal(i) / Decimal(2**75)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 75)
plot_pq(*pqzo, 3 * 75 + 32*75)
plot_pq_range(*pqzo, 50, 100)
plot_pq_range(*pqzo, 2500, 3 * 75 + 32*75)
plot_pq_range(*pqzo, 60, 75)
plot_pq_range(*pqzo, 2550, 2565)
In [7]:
getcontext().prec = 10000
epsilon = Decimal(1) / Decimal(2**100)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100)
plot_pq(*pqzo, 3 * 100 + 32*100)
plot_pq_range(*pqzo, 85, 100)
plot_pq_range(*pqzo, 93, 96)
In [8]:
getcontext().prec = 10000
epsilon = Decimal(1) / Decimal(2**100)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100 + 64*100)
In [4]:
getcontext().prec = 20000
epsilon = Decimal(1) / Decimal(2**50)
c = epsilon
a = epsilon
b = 1 - a - c
pqzo = make_pqzo(a, b, c)
plot_pq(*pqzo, 3 * 100, basic=False)
plot_pq(*pqzo, 6 * 100, basic=False)
plot_pq(*pqzo, 9 * 100, basic=False)
# Split up to prevent running out of ram
In [6]:
plot_pq(*pqzo, 20 * 100, basic=False)
# plot_pq(*pqzo, 40 * 100, basic=False)
# plot_pq(*pqzo, 60 * 100, basic=False)
# plot_pq(*pqzo, 80 * 100, basic=False)
# plot_pq(*pqzo, 100 * 100, basic=False)
# plot_pq(*pqzo, 120 * 100, basic=False)
# plot_pq(*pqzo, 140 * 100, basic=False)
# plot_pq(*pqzo, 160 * 100, basic=False)
# plot_pq(*pqzo, 180 * 100, basic=False)
# plot_pq(*pqzo, 200 * 100, basic=False)
# # plot_pq(*pqzo, 400 * 100, basic=False)
# plot_pq(*pqzo, 800 * 100, basic=False)
# # plot_pq(*pqzo, 1000 * 100, basic=False)
# # plot_pq(*pqzo, 1500 * 100, basic=False)
# # plot_pq(*pqzo, 2000 * 100, basic=False)

---------------------------------------------------------------------------
Overflow                                  Traceback (most recent call last)
<ipython-input-6-b2a85c60657d> in <module>
----> 1 plot_pq(*pqzo, 20 * 100, basic=False)
      2 # plot_pq(*pqzo, 40 * 100, basic=False)
      3 # plot_pq(*pqzo, 60 * 100, basic=False)
      4 # plot_pq(*pqzo, 80 * 100, basic=False)
      5 # plot_pq(*pqzo, 100 * 100, basic=False)

<ipython-input-3-3d26527cd5ca> in plot_pq(p, q, z, o, k, **kwargs)
      1 def plot_pq(p, q, z, o, k, **kwargs):
----> 2     plot_pq_range(p, q, z, o, 0, k, **kwargs)
      3 
      4 def plot_pq_range(p, q, z, o, k_low, k_high, simplex=True, simplex2=True, basic=True, log_scale=False):
      5     if simplex:

<ipython-input-3-3d26527cd5ca> in plot_pq_range(p, q, z, o, k_low, k_high, simplex, simplex2, basic, log_scale)
     11             y_axis = []
     12             for k in range(k_low if fun != o else max(k_low, 1), k_high):
---> 13                 point = np.array([fun(1, k), fun(2, k), fun(3, k)])
     14                 x_axis.append(point @ X)
     15                 y_axis.append(point @ Y)

<ipython-input-2-a40563d1c2bd> in p(i, k)
      6         if k == 0:
      7             return [p1, p2, p3][i - 1]
----> 8         return p(i, k - 1) * (p(i, k - 1) + 2 * p(i % 3 + 1, k - 1))
      9 
     10     @cache

Overflow: [<class 'decimal.Overflow'>]
In [7]:
def plot_trajectory(abc, k):
    a, b, c = (Decimal(x) for x in abc.split())
    pqzo = make_pqzo(a, b, c)
    plot_pq(*pqzo, k, basic=False)
In [8]:
plot_trajectory(".3 .4 .3", 100)
In [9]:
plot_trajectory(".1 .8 .1", 100)
In [10]:
import random
random.seed(0)
for i in range(5):
    a = random.randrange(1000)
    b = random.randrange(1000 - a)
    print(f".{a:03} .{b:03} .{(1000-a-b):03}")
    plot_trajectory(f".{a:03} .{b:03} .{(1000-a-b):03}", 100)

.864 .098 .038
.776 .107 .117
.041 .265 .694
.988 .008 .004
.497 .207 .296